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Abstract 


This  paper  studies  how  sampling  variation  in  component 
reliability  estimates  affects  the  computation  of  system 
reliability  that  uses  these  estimates  as  input.  Results  show 
that  relative  bias  in  system  reliability  grows  quadr a t i c al  1  y  with 
the  number  of  components  for  which  each  component  reliability 
estimate  is  used,  whereas  the  corresponding  coefficient  of 
variation  grows  linearly  with  this  number  of  components.  If 
these  components  are  in  parallel  they  lead  to  an  understatement 
of  system  reliability.  In  series,  they  lead  to  an  overstatement. 
The  paper  describes  resampling  schemes  that  eliminate  bias 
without  increasing  the  dominant  variance  term. 
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Introduction 

Every  computation  of  system  reliability  relies  on  the 
availability  of  numerical  values  for  the  reliabilities  of 
components  from  which  the  system  is  constructed.  If  these 
numerical  values  were  exact,  then  a  direct  computation  of  system 
reliability  would  at  most  suffer  from  numerical  roundoff  error. 
Since  the  numerical  values  of  component  reliabilities  rarely  are 
known  with  exactness,  a  system  reliability  computation 
customarily  employs  estimates  of  these  quantities  derived  from 
test  data.  This  substitution  introduces  an  additional  source  of 
error  attributable  to  the  sampling  variation  inherent  in  the 
component  reliability  estimates.  As  the  present  paper  shows, 
neglecting  this  source  of  error  can  produce  a  misleading  system 
reliability. 

This  error  manifests  itself  in  bias  and  variance.  For  a 
system  composed  of  several  types  of  components  where  the  system 
reliability  computation  uses  a  common  component  reliability 
estimate  as  input  for  all  components  of  the  same  type,  the 
relative  bias  in  system  reliability  increases  qua  dr  a t i ca 1 1 y  with 
eacn  of  the  numbers  of  components  of  each  type,  whereas  the 
corresponding  coefficient  of  variation  grows  linearly  with  these 
numbers.  For  components  of  the  same  type  in  parallel,  this 
system  reliability  computation  understates  true  reliability.  For 
components  of  the  same  type  in  series,  the  computation  overstates 
reliability. 

These  results  imply  that  for  given  component  reliability 
estimates  system  reliability  computations  for  two  different 
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syatems  composed  of  exactly  the  same  number  of  components  of  each 
type  can  have  substantially  different  statistical  error 
characteristics.  While  no  method  exists  for  reducing  the 
variance  of  the  system  reliability  base  on  component  reliability 
estimates  of  fixed  sample  sizes,  r esamp 1 i ng  schemes  do  allow  one 
to  eliminate  bias  without  increasing  the  dominant  variance  term. 

Section  1  introduces  the  notation  for  characterizing  a 
system  as  a  network.  Section  2  gives  the  conventional  estimator 
for  system  reliability  and  describes  how  one  can  use  a  confidence 
interval  to  assess  its  statistical  accuracy.  Section  3  shows  how 
parallel  and  series  systems  affect  statistical  error  and  Section 
4  extends  the  results  to  more  general  systems.  Section  5 
describes  two  resampling  plans  that  eliminate  bias  while 
preserving  the  dominant  variance  term.  Section  6  gives  the 
conclusions  of  the  study. 

1  .  System  Char ac t er  i  s t  i  c s 

As  a  basis  for  studying  error,  consider  tne  network  G  =  (V,E) 
with  node  set  V  and  arc  set  E.  For  convenience  of  exposition, 
assume  that  nodes  represent  components  that  function  perfectly 
and  that  arcs  represent  components  that  fail  randomly  and 
independently.  Hereafter,  we  treat  the  word  component  as 
synonymous  with  arc.  To  characterize  G  more  completely,  we 


define: 


r  =  number  of  distinct  types  of  components 
p^  =  probability  that  a  component  of  type  i  functions 

i  =  1 . r 

p  =  (  P1 . Pr  ) 

E i  =  set  of  arcs  that  use  components  of  type  i 

r 

(E  OE  -0  i-j  ,  E  -  u  E  ) 

J  i-1 

k  =  1-^1  number  ofl  components  of  type  i 

k  »  (  k  . k  ) 

_  1  r 

e..  =  jth  arc  in  E. 

ij 

x.  .  =  1  if  arc  e.  .  functions,  =  0  otherwise 
i  J  i  J 

vi 

x.  -  I  x  *  number  of  arcs  of  type  i  that  function 


x  -  (  xn  ’  *  *  ‘  ’  x!  k1  s  X2 1  *  '  "  •  *  x2k2  : 
X  -  set  of  all  arc  states  x 


,  x  .  ) 

r  k 

r 


r  x  . 

P(x,k,p)  -  IT  p  1  (1-p  ) 
-  -  -  i-1 


kl-xl 


=  probability  mass  function  of  states  in  X 

<P  ( x )  =  1  if  the  system  functions,  =  0  otherwise 

g(p)  =  I  4>  (  x )  P(x,k,p)  =  probability  that  the  system  functions 
xeX  ~  '  ~ 

We  also  assume  that  G  describes  a  coherent  system.  A  system  of  com¬ 
ponents  is  coherent  if  its  structure  { 4>  (  x )  }  is  nondecreasing  and 
each  component  is  relevant.  See  Barlow  and  Proschan  (1981,  p.  6). 

The  system  reliability  g(p)  can  have  diverse  i n t er pr e t a t i ons . 

For  example,  let  T  denote  a  subset  of  V  and  let 

4>  ( x )  =  1  if  all  nodes  in  T  are  connected  when  arc 


state  x  occurs 


0  otherwise. 


Then  g(p)  denotes  the  probability  that  all  nodes  in  T  are 


r 


connected.  If  T  «  js,t},  this  is  called  the  s-t  connectedness 
problem.  If  T  =  V,  it  is  called  the  all  terminal  connectedness 
problem  . 

Reliability  in  flow  problems  can  also  be  characterized. 
Suppose  that  G  is  a  directed  acyclic  flow  network  with  source 


node 

s  and 

term! nal 

node  t 

Let 

pi 

=  pr 

(arc  j 

has  flow  capacity 

b  .  )  b . >0 

l  l 

1  - 

pi 

=  pr 

(arc  j 

has  zero  flow  capacity) 

Xij 

=  1 

if  arc 

j  in  Ej  has  flow 

capacity  b^ 

=  0 

if  arc 

flow  capacity  is 

zero 

and 

let 

<P  ( x  ) 

=  ip  (  X 

,  z ) 

=  1 

if  the 

maximal  s-t  flow 

exceeds  a  3oecified 

demand 

z  when  state  x  occurs 

=  0 

otherwise. 

Then 

g(p)  = 

g  ( P 

,  Z  ) 

denotes 

the  probability  that  the  maximal  s-t 

L 


i 


'A 


A 


t 


flow  in  G  exceeds  z. 

Although  the  exact  computation  of  g(p)  for  these  examples 
belongs  to  the  NP-hard  class  of  problems  (Valiant  1979,  Ball  and 
Provan  1983,  Provan  1986),  we  assume  that  for  a  particular 
network  instance  of  interest,  one  can  indeed  effect  the  exact 
computation  i_f  p  is  known.  If  an  exact  computation  proves  infeas¬ 
ible  and  one  resorts  to  the  Monte  Carlo  method,  then  one  needs  to 


51 
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perform  a  more  elaborate  analysis  to  determine  how  statistical 
variation  in  the  component  reliability  estimates  interacts  with 
the  sampling  variation  that  the  Monte  Carlo  method  induces. 


2 .  Component  Reliability  Estimates 

In  practice,  p  is  not  known  exactly,  but  can  be  estimated 
from  test  data.  Suppose  one  tests  components  of  type  i  for 
i  - 1  ,  .  .  .  ,  r  .  Each  test  begins  with  a  new  component  functioning. 
Let  Z  j  j  denote  the  outcome  of  the  jth  test  of  component  of  type  i 
where  Z^j  »  1  if  the  component  functions  at  the  end  of  the  test 
period  and  Z^j  =  0  if  the  component  fails  prior  to  the  end  of  the 
test  period.  Presumably  each  component  of  type  i  is  tested  under 
identical  conditions  that  resemble  the  system  environment.  Then 


one  has  the  data  vectors  Z^  »  {z^.-.-.Z 


for  i  ='  1  ,  .  .  .  ,  r  where 


the  elements  of  Z^  are  independent  and  identically  distributed 
with  p^  =  EZ^j  j=1 . n^,  are  independent,  and 


P.  =  nT1  I1  Z.  . 
1  1  J-1  1J 


gives  the  maximum  likelihood  estimator  of  p.  with 


Ep  .  =  p  . 

l  i 


var  p.  =  pi  ( 1 -p .  ) /n  t 


E  (  p  .  -  p  .  )  m  -  o(i/n  J— (m  +  1)/2J 

I  X  1 


)  m  =3  ,  **  ,  .  . 


as  n  . 


ysnr*  \m  w  wr. 


LtJMJ.U 


where  0 ( y )  as  y  +  L  denotes  a  function  f  such  that  lim  | f ( y ) j / y  is 

y+L 


bounded.  Observe  from  (2)  that  p.  is  an  unbiased  estimator  of  p.. 

1  l 


Let  p  =  ( p 1 . p^ )  .  Then  it  is  not  unusual  to  estimate  g(p) 


by  g(p).  Although  other  methods  exist  for  using  test  data  to 


estimate  component  reliabilities,  the  appeal  of  the  method  that  we 


adopt  here  arises  from  the  well-understood  sampling  properties  of 


p,  enabling  us  to  concentrate  on  the  statistical  variation  in  g(p) 


that  substitution  of  p  for  p  in  g(p)  induces.  As  Gaver  and  Hoel 


(1970)  show,  other  methods  can  lead  to  bias  in  component 


reliability  estimates,  which  would  force  us  to  conduct  a  more 


complicated  analysis  to  get  at  the  sampling  properties  in  the 


system  reliability  estimate. 


As  Sections  3  and  4  make  clear,  g(p)  generally  either  under¬ 


states  or  overstates  g(p)  with  regard  to  expectation.  However, 


for  the  moment,  we  describe  how  one  can  globally  assess  the 


statistical  accuracy  of  g(p)  based  on  confidence  intervals 


computed  for  p^,...,pr. 


Let  Z 


n  . 

I  Z.  ..  For  each  p.  we  seek  a  100*(1-a)  confidence 
1  j=1  1J 


interval  [p*(Z.,n.),  p**(Z.,n.)] 
ill  i  ii 


pr[p*(Z.,n.)  <  pi  S  p**(Z.,n.)]  > 


1  -  a. 


F  (  m  ,  q  ) 


l  (?)  Q  1  ( 1 - q  ) 1 

i  =0 


0 < q < 1  ;  j  =0 , 1  ,  .  .  .  , m  ;  m= 1  , 2 , 


"viva’ll 


F  (n,q_)  =  1  -  a/2  for  i=1,...,r 

Z  *— 

for  p  (n,z)  and  (n,z),  respectively,  and  achieve  a  confidence 
coefficient  of  at  least  1-a.  We  call  the  result  a  binomial 
interval  . 

As  nj  increases,  exact  solution  becomes  difficult  because  of 


numerical  error.  Then  one  has  the  well  known  result 

1 im  pr { 

n  .  ■»“ 

1 

where 

c(a)  =  {y:  (2tt)  ^  e  W  ^  2 dw  =  1-a/2 }, 

—  00 

and  in  principle  one  can  solve  the  corresponding  quadratic  form 
p 2 [ 1  +c 2 ( a ) / n .  ]  -  p . [ 2p ^ +c 2 ( a ) / n .  ]  +  p2  =  0 

for  p*(n.,Z.)  <  p**(n.,Z.).  The  resulting  confidence  interval 

ill  i  i  i 

an  error  of  approximation  which  decreases  as  increases. 
ever,  the  rate  of  convergence  is  nonuniform,  being  most  rapid 
Pi  *  */2  and  least  rapid  for  p ;  close  to  zero  and  unity. 


I  P  i  -  p  i  I  , 

— — - - - jt  £  c(a)}  =  1  -  a, 

[p  ( 1 -p  )/n  ]  2 


(6) 

has 

How 

for 


This  non- 


uniform  convergence  limits  the  appeal  of  this  confidence  interval 


in  practice. 

A  third  approach  uses  Chebyshev's  inequality  so  that 
p * ( n i , Z i )  £  p**(ni,Zi)  are  again  solutions  of  (6),  but  with  1 /c/2 
replacing  c(a).  Although  this  confidence  statement  holds  for 
every  n ^ ,  the  interval  width  can  be  wide.  A  fourth  approach 
based  on  the  probability  inequality  (Okamoto  1955,  Hoeffding 
1  963  ) 

p.+e  l-p.-en. 

pr(Zi~pi>e)  <  lCp./(p.+e)]  1  [ ( 1  - p i ) / ( 1  - p i - e  )  ]  1  }  1 

0<e  < 1  -  p . 

l 

produces  tighter  intervals  for  small  a.  For  n ^ Sin ( a/ 2 ) / In  max(p. ,1- 
(pj^n^.Z^),  p**(n^,Z^))  covers  p^  with  probability  >  1  -  a  where 
p * ( n i , Z i )  <  p**(n^,Z^)  are  now  the  solutions  to  the  equation 

Pt  1 n ( 0 / p .  )  +  (1-p.)  ln[(1-0)/(1-p.)  ]  =  J  1 n ( a/ 2  )  .  (7) 

See  Fishman  (1986).  The  ratio  ln(a/2)/ln  max(p^.l-p^)  provides 
an  indication  of  whether  or  not  n j  exceeds  the  required  lower 
bound  . 

Although  the  resulting  interval  leads  to  a  confidence 
interval  of  greater  width  than  the  binomial  and  normal  intervals 
do,  it  is  considerably  easier  to  compute  than  the  binomial 
interval  is  for  moderate  and  large  n  i  and  induces  no  error  of 


approximation  as 


the  normal  interval  does. 


Therefore,  we 


recommend  the  computation  of  a  binomial  interval  from  (5)  when. 


possible  and  an  interval  based  on  (7)  otherwise,  provided  that 
n  i  >  ln(a/2)/ln  max  (  pj_  ,  1 -p  ^  )  . 

Since  . are  independent,  one  has 

r 

pr  jpe  JI  [p*(n  .  ,  Z.  )  ,p**(n  .  ,  Z  .  )  ]  }  >  B 
i » 1  1  1 

where  6  =  (  1  -  a  )  r  .  To  achieve  a  confidence  level  B,  one  chooses 
a  =  1  -  B1/r  for  each  interval.  Since  the  system  is  coherent, 

3g(p)/3pi  £  0  for  i  =  1,...r.  Therefore 

pr[g(p)<g(p)<g(p  ) ]  £  B  (8) 

u  v  u  M.  it  li  X  X 

where  p  =(p1(n1>Z1) . p*(nr,Zr>)  and  p  =  (p]  ( n 1  , Z 1  pp 

Since  p  * ( n  ^ , Z  ^ )  S  p  i  S  p  *  ( n  ^  ,  Z  ^  )  for  i  =  1 . r  with  probability 

one  and  since  3g(p)/3p.  >  0  for  i=1,...,r,  one  has 

*  -  #  # 

g(p  )  2  g(p)  Z  g(p  )  with  probability  1, 


a  result  which  provides  a  convenient  way  of  assessing  the  extent 
of  sampling  variation  in  g(p).  With  p*  and  p**  in  hand  for 
specified  B,  one  can,  for  a  specified  system  G,  compute  g(p*)  and 

g(p**)  in  two  reliability  evaluations  and  determine  whether  or  not 

*  *  * 

the  interval  width  g(p  )  -  g(p  )  is  sufficiently  small  for  the 
purposes  of  reliability  analysis.  As  Sections  3  and  A  show, 
there  is  good  reason  to  believe  that  this  interval  grows 


substantially  as  the  size  of  the  system  G,  constructed  from 


components  of  types  1,...,r,  grows . 

3 .  Parallel  and  Series  Systems 

We  use  the  s-t  connectedness  problem  to  illustrate  the 

a 

potential  seriousness  of  errors  in  the  estimate  g(p). 

Theorem  1  .  Let  G  denote  a  network  of  k-|  arcs  of  type  1  in 

parallel  with  source  node  s  and  terminal  node  t  so  that 

k1 

g(p)=1-(1~P1)  (9) 

gives  the  probability  that  s  and  t  are  connected.  Let 

denote  0-1  test  data  on  ni  components  of  this  type, 

1  1  n  1  A  ^  k  1 

let  Pi  =  -  l  Z  ,  p  =  p  and  g(p)  »  1  -  (1-p  )  .  Also 

1  j  =  1  1 J  -  1 

Chebyshev's  inequality  gives 

pr{|p1-p1|<B[p1(1-p1)/n1]/z}  >  1  -  1  /  6 2  (10) 

8  >0  . 

Based  on  (10),  the  minimal  sample  size  required  to  achieve 

pr  1  | g ( P  )~g( P  )  | < e [ 1 -g ( p  )  ]  }  >  1  -  -U  e>0  (11a) 


(11b) 


n .  >  62p,/(1-p1)[(1+e) 


1  /k 


’-1  ]2 


and 


*  2 

1  im  n1 / k 1 

k  -*•« 


>  B2P1/C(1~P1 )Cln(1 +e)]2 


(11c) 


Proof .  Substitution  into  (11a)  gives 

1/k 

pr { | g ( P ) -g ( P ) | < e [ 1  - g ( P ) ] }  =  pr{[1-(1+e)  ]  ( 1  - p 1  ) <p 1 - p 1 

1  /ki  i 

<  [  1  -  ( 1  -  e  )  Kl-p,  )  }  >1  -  J2- 

Then  Chebyshev's  inequality  (10)  gives 

2  ^  ^ k  1  2  ^  ^ k  1  2 
n}  >  B  ^p 1 / ( 1  - p 1 )  mi n { [ 1  - ( 1 +e  )  Cl-(l-e)  ]  }. 

1  /k1  1 /k1  * 

Since  (1+e)  -1  <  1  -  (1-e)  ,  n.  in  (11b)  follows.  Expression 

1/k 

(11c)  follows  by  applying  L'Hopital's  rule  to  (1/k^)/[(1 +e)  -1] 

Note  that  (11a)  is  an  attempt  to  control  the  relative  error 
on  the  system  failure  probability  1-g(p).  Expression  (11c) 
immediately  makes  apparent  that  the  sample  size  n*  needed  to  keep 
this  relative  error  at  e  grows  quadratically  with  k i ,  the  number 
of  arcs  in  parallel.  Theorem  2  provides  insight  into  the  source 


of  the  potential  error. 


Theorem  2. 


For  k1  arcs  of  type  1  in  parallel  with  g(p1) 

a  A  k « 

1  -  (1-p^  and  g(p)  =  1  -  ( 1  - p  1  )  , 


lim  Eg ( p )  -  1  -  ( 1  - p  ) 
k^® 


k1_1 

lim  n1 E[g(p)-g(p)  ]  =  - k 1 ( k 1 -  1 ) p 1 ( 1  - p 1  )  /2 


and 


*  2 

lim  n 1  var  g(p)  =  k1  p, (1-p] ) 


Proof .  Since 

Eg ( p  )  =  1 


(1 


Pi) 


n .  - 1  k  i. 

-  'l  <J/n.)  '  (,') 
j  =0  J 


i  nrJ 

p, ( i -Pi ) 


(12)  follows  immediately.  Let  A  =  p  i  -  p  -j  and  observe  that 

„  k  k.k.  k  —  j  j  j 

g(p)  -  1  -  ( 1  - p i - A )  =  1  -  l  [  ]  ( 1  - p  i  )  A  (-1) 

J-0  J 

k i  k1_^  ^ i “ 2  * 

*  1  -  (  1  ~  P  i  -  A  )  +k  1  ( 1  -  p  1  )  A-k  1  (  k  i  -  1  )  (  1  -p1  )  A 

Since  EA  =  0,  EA2  =  p1(1-p1)/n1  and  EAm  =  0  (  1  /  n  1l-m  + 1  ^  7  ^  )m 

k  1  -  1  2 

Et-gCp)-g(p)  ]  =  -  k  i  (  k  i  -  1  )  p  1  (  1  -  p  1  )  /  n  i +0  (  1  /  n  j  )  as  n,-* 

and  (13)  follows.  An  analogous  development  gives  (14). 


(12) 

(13) 

(14) 


/  2  +. 

£  3 , 


Theorem  2  raises  several  concerns.  The  quantity  g(p)  under- 


E 


states  the  true  reliability  g(p).  Moreover,  the  relative  bias 

E[g(p)-g(p)  ] 

i  -  g<p)  -  kt(kr° 

makes  clear  that  the  dominant  term  in  the  relative  understatement 
increases  quadr a t i cal 1 y  in  kj.  If  the  objective  is  to  design  a 
parallel  system  based  on  component  of  type  1  with  a  specified 
level  of  system  reliability,  then  g(p)  encourages  one  to  add  more 
components  in  parallel  than  may  truly  be  required. 

Observe  that  the  coefficient  of  variation 


Y  ( k  ,  p  ,  n  )  »  [var  g  (p  )  ]/z/[  1  -g  ( p  )  ]  -  k  1  [p  1  /  ( 1  -  p  1  )  n  1  Y1 


reveals  linear  growth  in  k^.  As  a  result,  a  network  with  Jk^ 
components  of  type  1  in  parallel  would  lead  to  a  coefficient  of 
variation  J  times  larger  than  a  network  G  ]  with  just  k -|  com¬ 
ponents  in  parallel. 


An  analogous  development  for  k^  >  1  components  of  type  1  in 

~k  1 

series  gives  a  sample  system  reliability  g(p)  *  P1  that  over¬ 


1 


states  the  system  reliability  g(p)  =  P1  .  Again,  relative  bias  is 


proportional  to  k^  and  the  coefficient  of  variation  is  proportional 


to  k.j.  More  generally,  consider  a  set  of  r  subsystems  in 
series  where  subsystem  i  is  composed  of  kj  >  I  components  of 


■  J*  •  •  •  •  -  •  •  ■  -  '  »  ‘  *  JI  V  V«  J  ^  V  «  *  «  -N  -N  %  .  »  .  »  \  • .  V  ■ 


y  .■ 


type  i  in  parallel  i  =  1  . r.  Here  system  reliability  is 

r  k  ^  r  k 

g(p)  -  II  [  1  -  (  1  -  p  .  )  1]  and  clearly  the  quantity  g(p)  =  n  [  1  -  (  1  -  p  1  )  ] 

i-1  i=l  1 

understates  it.  Conversely,  a  set  of  r  subsystems  in  parallel 

where  subsystem  i  has  k^  >  1  components  in  series  has  reliability 

r  k .  ^  r  ^  k . 

g(p)  -  1  -  II  (1-p.)  and  the  quantity  g(p)  =  1  -  n  (1-p.1  )  over- 

i  =  1  1  ~  i=1  1 

states  it. 


i 


i 

I 


» 

I 


i 

r 

i 

! 

I 

« 

! 

t 

* 

* 


! 


i 

5 

) 
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4 .  More  General  Systems 

Results  for  more  general  systems  reveal  how  potential  errors 
grow  with  the  number  of  types  of  components  r  as  well  as  with  the 
number  of  components  of  each  type. 

■  Theorem  3  •  Consider  a  system  composed  of  . . . kr  components  of 

types  1 ,  .  .  .  ,  r  respectively.  Then 

Eg (p )  =  g(p)  +  w(k,p,n)  +  R1  (15) 


and 

var  g(p)  =  v (k , p ,n )  +  R2>  (1 6) 


where 


r  k(k.-1)p-2x(k-1)p.+x.(x-1) 

w(k,p,n)  *  l  {  l  <j>(x)P(x,k,p)[— — - - — j - - — — — - — - 

i=1  xeX  -  p  (1-p.  )c 


Pi(1'Pi} 
2n . 


(17) 


v 


(18) 


r  x.-k.p.  _  p. (1-p. ) 

v(k,p,n)  =1(1  4>(x)  P(  x ,  k  ,  p )  [-1,  1  1  J}  — — 

-  i=1  xeX  -  pi  i 


and 


i 

l  0(1/n  n  ) 
i,j-1  J 


as  min  n  -*■« 
1£i<r  1 


R?  -  l  0( 1 /n. n . ) , 
i  ,  j  =1  1  J 


as  min  n .  -*•” 
1 <i<r  1 


Proof  .  Let  A  ^  =  p  ^  -  p^  and  observe  that 


g(p)  -  l  4>(x)  P(x,k,p) 
xeX  "  - 


=  I  *<x)  n  [  r 

xeX  ~  i=1  j=0 


•x.  -m 


k  -x.  x.  k  -x.  x.-i  ...  . 

I  1(i1)(  _  M  p,1  (i-P,)  1  1  H)V+m] 

m=0  J  1  1  1 


=  I  <t>(x) 
xeX 


r  x.  k.-x.  x.  k.-x.  x. 

n  [p/h-pj  1  1  vi  '  ' 


lx  xl  i(7)(ki"Xi) 


i  =1 


j  =0  m=0 


J 


(-1 )mA^ +m 
_ I 

p^ ( 1 -p . )m 

l  l 


]• 


Expressions  (15)  and  (16)  following  from  substitution  of  (2),  (3) 

and  ( 4 )  for  EA^  m  for  j  ,  m =0 , 1 . k  - x  and  the  observation  that 

that  EA.A.*=  0  for  i*i’. 
l  i 

2  2 

In  addition  to  the  proportionality  to  k . k  ,  observe  that 

l  r 

the  number  of  terms  in  w(k,p,n)  and  v(k,p,n)  increases  linearly 
with  r,  the  number  of  types  of  components.  This  increase  would 
become  quadratic  if  the  data  vectors  Z  ^  ,  .  .  .  ,  Z  were  positively 


correlated. 


alternative  method  of  using  the  current  data  more  effectively. 

Unfortunately,  variance  reduction  is  not  possible.  Since  p 

is  the  maximum  likelihood  estimator  of  p,  g(p)  is  the  maximum 

likelihood  estimator  of  g(p)  and  v(k,p,n)  corresponds  to  the 

Cramer-Rao  lower  bound  on  variance  for  fixed  k  and  p  as 

n  1  ,  .  .  .  ,  n  -»« .  That  is,  no  alternative  estimator  of  g(p)  based  on 

Z, can  achieve  an  asymptotic  variance  smaller  than  v(k,p,n) 
l  r  ~ 

in  ( 1 8  )  . 

The  potential  for  bias  removal  is  more  promising.  Recall 
that  positive  bias  can  lead  to  a  more  frequent  failure  pattern  in 
practice  than  the  computed  reliability  implies.  A  negative  bias 
can  lead  to  a  costly  enhancement  of  the  system  to  mitigate  the 
apparent,  but  not  real,  reliability  deficit  that  the  reliability 
computation  suggests.  This  section  describes  a  method  of  removing 
this  statistical  bias  while  preserving  the  asymptotic  variance  at 
its  minimum  v(k,p,n).  The  method  uses  a  data  resampling  plan  to 
produce  an  unbiased  estimate  of  system  reliability  in  time  per 
trial  that  grows  considerably  more  slowly  than  the  time  required 
to  compute  g(p). 

Recall  that  data  vectors  Z,  ,  .  . .  ,  Z  which  were  used  to  esti- 
mate  pj,...,pr  and  assume  that  n j  >  ki  for  i  =  1,...,r.  Algorithm 
A  describes  a  procedure  that  on  each  trial  (step  2)  randomly 
samples  (without  replacement)  and  assigns  an  element  of  Z.  to 
each  component  of  type  i.  Let  X  denote  the  resulting  arc  state 
vector  of  zeros  and  ones.  Given  this  assignment,  the  system 
either  functions  ( $ ( X  )  = 1  )  or  fails  ($(X)= 0).  Then  hK  (step  3) 
is  our  refined  measure  of  system  reliability. 


Algorithm  A 


Purpose:  To  compute  an  unbiased  estimate  h „  of  system 

1  J\ 

reliability  g ( p ) . 

Input :  Network  G  -  (V,E),  where  E  *  { e  1  1  ,  .  .  .  , e 1 k 

er  1 . e^k  1»  sample  data  Z  ^  =*  (Z^^F*..FZ^^ 

r 

i  =  1 . .  and  desired  number  of  trials  K 

Output:  h„. 

-  I\ 

Nomenclature:  X  =  (X . ,X,,  ;...,X  X  ,  ). 

-  _  11  lk,  r l  r  k 

1  r 

Method  : 

1  .  Set  S+- 0  . 

2.  On  each  of  K  trials: 

a.  For  i  *  1 , . . . , r  : 

1 1  i  •  •  •  i  n  j  |  • 

For  j  »  1,...,kj:  sample  e  from  ;  remove 

e  from  W.;  set  X  .  .  Z  ,  . 

-i  i  J  i  e 

b.  Determine  <p  (  X  )  ;  set  S<-  S  +  <|>  (  X  )  . 

3.  Compute  reliability 

h  S/K. 

K 


End  of  procedure 


Theorem  4 . 


For  hK  as  computed  in  Algorithm  A, 


a  . 


b  . 


c  . 


EhK  =  g(p) 


var  hK  =  g ( p ) [  1 -g ( p ) ]/K  +  [v(k,p,n)+  l 

i  = 
r 

lim  var  h  =  v(k,p,n)  +  I  0(1/n.n.) 

K  —  K  -  i,j=1  1  J 


r 

l  0(1/n.n.)](K-1)/K 
J  =  1  1  J 

as  min  n  .  . 

I <i  <r  1 


Proof  of  a.  Observe  that 


r 

<p(X )  =  l  <p(  x )  n 

xeX  ~  i = 1 


k  . 

n1 

j=i 


(i-x.  . ) 

i  j 


x  .  . 
1 J 


where  the  X^j's  are  sampled  in  step  2a.  Since  sampling  occurs 
without  replacement  on  each  trial 


E  ^[X^O-X..)1"^ 
J-1  1J  1J 


]  = 


11 1  ECX^O-X..)1"5^ 

J-1  1J  1J 


] 


x  .  ,  .  .  k  .  -  x  .  . 

=  p  .  l  (  1  -  p  )  l  l 

l 


Theref  ore , 


E<D(X)  =  g(p) 


and  consequently  h  is  unbiased. 

l\ 


Also, 


since  4>  (  X  ) 


4>  (  X  ) 


var  <f»  (  X  )  =  g(p)[1-g(p)  ]. 


-2 


Proof  of  b.  Let 


Q  -  ( q i i  » • • • »  q i k i  ;  •  • •  :  Qr i . Qrk^  ) 


and  redefine  the  reliability  function  as 


r  k  x  i  —  x 

h(q)  =  l  4>(x)  n  n1  q.^(1-q.  )  i'i 
xeX  ~  i  =1  j=1  J  J 


Now  write  the  Taylor  expansion 


h ( q )  =  h(p  . p;...;p,...,p)  +  l  I 1 3  h / 3  q 

-  '  '  ”  i  .■  1  •  1  1  J 


j,  .Ij3h/Jdul  ,-p  <lij"pi)  *  a 

1=1  J  = 1  1J  1 


where  R  denotes  the  remainder  composed  of  higher-order  cross- 

(  V  )  (  z  ) 

derivatives.  Let  X.  and  X. .  denote  the  assignments  for  arc 

i  J  iJ 

ei  .  on  trials  v  and  w  respectively.  Then  for  all  j i  , j 2  =  1 . k 


ECU^-PjHX^’-p^]  -  Pjd-Pj) 


if  y  =  z  and  j  -j  =  j  2 


Then 


=  p .  ( 1  - p .  )  / n .  if  y  *  z 

r  l  l  l 


otherwise . 


Y(y)  fy(y)  Y(y)  .  Y(y)  Y(y) 

1  1  1  .  Ik.  ’  ■  *  *  ’  rl  .  rk 

1  r 

4<y)  -  x(y)  - 

ij  ij  ij 


h  (  X  f  y  5  ) 


p  k  x  I  —  x 

l  $(x)P(x,k,p)  II  Hi[1+A^y^(2x..-1)/p  1  ^  (  1  -  p  .  )  1 J  ] 

x  e  X  '  "  “  "  i=1  j=1  1J  1J 


l.»  U*  .!*•  4#»J 


4**1 


-21  - 


Since 


8g/3p.  =  I  3h/8q.  . 

-  i  .  .  l  j  1  q  .  . =p  . 

J=1  ij  i 


one  has  for  y*z  . 


cov[h(X^y)),  h ( X  ^  Z  ^  )  ]  =  l  ( a  g / 3  p . ) 2  p  (l-p.)/n. 

i  =  1  iiii 


+  l  0 ( 1 /n . n . ) 
i  ,  j  =1  1  J 


as  min  n .  -*•“ 
1 <i <r  1 


Since 


v ( k  ,  p  ,  n  )  =  l  ( 9g/3p  ) 2  p  (1-p.)/n. 
"  "  “  i-1  1  11 


as  min  n  .  , 

iSiSr  1 


the  quantity 


*  K  =  R 

y-i 


var  h. 


var  h(X(y))/K  +  cov[h(X(y)),  h(X(z)  )  ]  ( K  —  1  )/K 


=  g(p)[1-g(p)]/K  +  [ v ( k , p  ,  n  )  +  l  0(1/n.n.)](K-1)/K 

-  -  -  i  ,  j=i  1  J 

as  min  n  .  ■*<*  , 

1 <i <r  1 

which  proves  part  b.  Part  c  follows  immediately. 

The  significance  of  Algorithm  A  is  now  apparent.  The 
resampling  scheme  produces  an  unbiased  estimate  of  g(p).  As  the 


number  of  trials  K  increases,  the  variance  of  h 


K  converges  to 


a  quantity  whose  dominant  term  is  the  Cramer-Rao  low 


bound  v ( k , p , n ) .  Moreover,  in  place  of  a  direct  calculation  of  the 
reliability  g(p),  one  computes  <f>  (  X  )  ,  in  step  2b,  K  times.  For  s-t 
connectedness,  a  depth-first  search  algorithm  computes  4>  ( X )  in 
0(max(|v|,|E|)  time.  See  Aho,  Hopcroft  and  U 1 1  mar.  (197  4).  If  G  is 
a  directed  acyclic  flow  network  with  random  binary  capacities  and 

<p  (  x  )  =  <(>  (  x  ,  z  )  =  1  if  maximal  s-t  flow  >  z 

=  0, 

one  can  determine  $ ( X )  in  0(|v|log|v|)  time  if  G  is  planar.  See 
Itai  and  Shiloach  (1979).  The  fastest  known  algorithm  for  a 
nonplanar  network  takes  0(|vj3).  see  Malhotra,  Kumar  and 
Maheshwari  (1978).  These  time  complexities  make  clear  that  the 
cost  of  resampling  per  tri*;  is  generally  incidental  relative  to 
the  cost  of  performing  the  exact  computation  of  g(p). 

A 

To  bring  var  hK  to  the  neighborhood  of  v(k,p,n)  one 
needs  to  make  K  sufficiently  large  to  make  g(p)[1-g(p)]  small 
relatively  to  v(k,p,n).  To  assess  when  this  occurs,  one  would 
need  to  observe  var  h  as  a  function  of  K.  This  quantity  is 

f\ 

unknown;  moreover,  it  is  not  possible  with  the  sampling  scheme 

of  Algorithm  A  to  compute  a  useful  estimate  of  var  hK. 

One  partial  solution  to  the  problem  partitions  the  data. 

Let  m i , . . . , mr  denote  integers  such  that  m  ^  >  k  ^  for  i=1,...,r,  and 

let  c =n i / m ^  =  .  .  .  =n r / m r  =  integer.  Then  Algorithm  B  describes  an 

alternative  scheme  that  involves  resampling  K*  times  from  each  of 

c  partitions  of  the  data  Z, . Z  .  Theorem  5  reveals  the  benefit 

- 1  _r 


of  this  method. 


Algorithm  B 


Purpose :  To  compute  an  unbiased  estimate  h  of  system  re¬ 

liability  g ( p )  and  an  unbiased  estimate  of  var  h. 

—  K. 

Input  :  Network  G  =  (V,E)  where  E  =  [e^  . e^  ;•••; 

e  1,...,erk  sample  data  Z.  =  (Z. ( . Z_n  ) 

r  ~  i 

i=1,...,r,  integers  c,m^f...tm  and  desired  number 

of  replications  per  partition  K*. 

Output:  h ,,  and  V  (  h  „  )  as  unbiased  estimates  of  g ( p )  and  var 

K  K  ~ 

Nomenclature :  X  =  ;...,X1,,..,Xrk  ). 

K 1  r 

Method  : 

1.  Set  K  «-  0  ;  For  y  =  1 . c:  set  Sy  «-  0. 

2.  On  each  of  K  *  trials: 

For  y  =  1  ,  .  .  .  , c  : 

For  i =1  ,  .  .  .  , r  : 

«-  {l,...,m^};  For  j=1,...,ki:  sample  e  from 

remove  e  from  H,  ;  set  X.  .  *•  Z.  ,  ,  s 

—  l  l  j  i , ( y- 1  )m  +e 

Determine  <(>  (  X  )  . 

Set  S  <-  S  +  (0  (  X  )  . 

y  y  _ 

3.  K  +■  K  +  c  K  * . 

4.  Compute  summary  statistics 


hK  <-  (  S1  +  .  .  .  +Sc  )  /K  . 


V<V 


c  (  c  -  1  )  y  =  1 


l  (cS  /K-h 
y  K 


End  of  procedure 


Theorem  5.  For  the  resampling  scheme  in  Algorithm  B, 


a .  Eh^  =  g(p  ) 

r 

b.  var  h  =  g ( p ) [ 1  - g ( p )  ]/K  +  v(k,p,n)  (K-c)/K  +  c  l  0(1/n.n.) 

K  -  -  -  i,j=1  1  J 

as  min  n  .  ->co 
1<i<r  1 

c.  EV ( h^ )  =  var  h^  . 

Proof.  Within  any  partition  y,  the  resampling  scheme  is 

identical  with  that  of  Algorithm  A  except  that  sampling  occurs 

from  Z.  ,  . . Z.  for  i=1,...,r.  Therefore,  for  each 

i , (y-1 )m.+1  i , ym . 

y  =  1 . c 

E(Sy/K*)  =  g(p)  , 

establishing  part  a.  Also 

p 

var ( S  /  K  )  =  g(p)[1-g(p)]/K  +  [  l  (  3  g  /  3  p  2  )  p  ( 1  - p  .  )  '  m 

y  ~  i  =  1  iiii 

p  *  * 

+  I  0  (  1  /  m  .  m  .  )  ]  (K  —  1  )  /  K  as  min  m.-*-°°  . 
i,j=1  1  J  1  <i  <r  1 

Since  Si , . . . ,SC  are  independent,  one  has 

var  h (  p  )  =  var(S  /K  )/c. 

K  1  y 

* 

Using  this  result,  together  with  m  =  n./c  i = 1 , . . . , r  and  K  =  cK  , 
gives  part  b.  Part  c  follows  by  taking  expectations. 


The  quantity  V ( h  )  provides  a  useful  estimate  of  var  h 

r\  \ 

which  one  can  use  sequentially  to  estimate  when  the  quantity 

g ( P ) [  1  ~g(  P ) ]/K  becomes  relatively  incidental  to  the  variance. 

That  is,  the  organization  of  Algorithm  B  enables  one  to  iterate  on 

step  2  to  generate  successive  estimates  ^cK*t  ^2cK*’’*'  and 

V(h  „#),  V(h_  „*),...  and  observe  the  extent  to  which  this 
cK  2c  K 

variance  measure  stablizes  as  a  function  of  K. 

The  one  drawback  of  Algorithm  B  as  compared  to  the  Algorithm 

A  arises  from  the  increased  relative  importance  of  the  higher 
r 

order  terms  £  0(1/n.n.).  These  are  scaled  by  c  in  Algorithm  B. 
i  .  j  =1  1  J 

As  the  sample  sizes  n] , . . . , nr  increase,  these  terms  diminish  in 
importance  in  each  case,  although  they  always  remain  c  times 


larger  in  Algorithm  B.  In  practice,  as  m -| 


mr  increase  c 


decreases,  reducing  the  significance  of  the  higher-order  terms. 
However,  a  smaller  c  means  a  less  statistically  reliable  estimate 


V(  hi/  )  of  var  hK  . 


6.  Conclusions 


In  general,  the  observations  made  in  this  paper  are  not 
encouraging  about  the  statistical  accuracy  of  a  system  reliabi¬ 
lity  computation  whose  input  consists  of  component  reliability 
estimates.  Although  no  alternative  system  reliability  estimator 
produces  a  smaller  asymptotic  variance,  the  resampling  schemes  of 
Section  5  do  provide  a  way  of  reducing  bias.  Based  on  the 
material  presented  here,  a  constructive  approach  to  system 
reliability  error  assessment  follows  these  steps: 

1.  Compute  component  reliability  estimates  p-|,...,pr. 

2.  Compute  1  00  x  (  1  -a  )  1  /'r  confidence  intervals  for  each 

component  reliability  p£  for  i=1 . r. 

3.  Compute  a  system  reliability  estimate  using  p^...^^  as 
input. 

4.  Compute  a  100*(1-a)  confidence  interval  for  system  reliabi¬ 
lity  using  the  confidence  intervals  for  . . . pr  in  step  2. 

5.  If  the  interval  width  for  system  reliability  is  within 
acceptable  bounds  at  coverage  level  1  -  a  ,  proceed  with  the 
study.  Otherwise: 

a.  One  may  improve  the  statistical  accuracy  of  the  point 
estimator  by  employing  the  resampling  schemes  in 
Section  6 

o  r 

b.  One  may  wish  to  collect  more  test  data  to  improve  the 
component  reliability  estimates  and  thereby  shorten 


the  interval. 
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